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Abstract 


The goal of this study was to experimentally and 
theoretically investigate the role of gravity-induced free 
convection upon the liquefaction dynamics of a cylindrical 
paraffin slab under normal gravity conditions. 

The experimental equipment consisted of a test cell, 
a fluid-loop heating system, and a multipoint recorder. 

The test chamber was annular in shape with an effective 
radius of 1.585 cm and a length of 5.08 cm. The heating 
chamber was a 1.906 cm diameter tube going through the center 
of the test chamber, and connected to the fluid loop heating 
system. 

By means of the heating chamber step changes in the 
inside wall temperatures were imposed on the test chamber. 
Temperature responses were measured by means of sixteen 
iron-constantan thermocouples mounted in the test chamber; 
hot wall temperature was also measured. N-octadecane , a 
high chain normal paraffin with an even number of carbon 
atoms, was used as the test material. All experimental runs 
were made with the longitudinal axis of the test cell in the 
vertical direction to insure that convection was not a func- 
tion of the angular axis of the cell. Ten melting runs were 
made (five sets of reproduced data) at various hot-wall 
temperatures. Also, two pure-conduction solidification runs 
were made to determine an experimental latent heat of fusion. 

The physical system was modeled using the energy equa- 
tion and the ideal-viscous flow velocity approximation. An 
Implicit Alternating-Direction (IAD) technique was used to 
approximate the partial differential equations in the com- 
puter solution. The IAD technique was used to eliminate the 
time step stability requirement inherent in an explicit fin- 
ite difference solution. However, because of numerical 
dispersion effects introduced by ignoring the second order 
partial in the Taylor's series expansion with respect to 
time, a time restriction was placed on the implicit solution 
in order to minimize the dispersion error. A latent heat of 
fusion of approximately 0.75 literature value was used in 
the theoretical solution. This value was determined by using 
a one-dimensional pure conduction model to predict the tem- 
perature results of the pure conduction solidification runs; 
the latent heat was varied until a match was obtained 
between data and theory. This new value was then used in 
the IAD program to model liquefaction. 



Experimental runs were made with the total temperature 
gradient of the liquid phase varying from a low of l6.5°K to 
a high of 27 . 8°K. There was very good agreement between data 
and theoretical prediction at the low temperature gradient 
levels. As the gradient became larger the deviation between 
theory and data became larger; all phase change times were 
predicted correctly, but final temperatures were not. 

The study has shown that with further refinements of 
the ideal-viscous flow model we should be able to predict 
temperature responses and phase-change times accurately; but 
the present model can only be used for theoretical predic- 
tion when liquid phase temperature gradients are small. 
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Introduction 


Phase' change thermal control techniques have received 
increasing attention (references 5, 6, and 7) in the last 
several years for spacecraft thermal design. Because of 
inherent advantages of simplicity and reliability a passive 
solid-liquid phase change material can be used in the walls 
of spacecraft as packaging around sensitive electronic 
equipment to absorb or release energy to maintain constant 
temperature of the electronic equipment. However, this 
system is limited by the heat rejection or absorption capac- 
ity of the material used. 

A previous study (lj has determined the property 
requirements of phase change materials in order that they 
be good thermal control devices. The material should be 
non-toxic, chemically-inert and stable, noncorrosive, have 
small density variations, and have a high latent heat of 
fusion. The material should also melt in 50- to 150°P 
range; n-paraffins with an even number of carbon atoms are 
the most widely used materials for this purpose. In this 
study, n-octadecane was used as the test material. 

An earlier study (2) at the Colorado School of Mines 
dealt with the problem of gravity-induced free convection in 
the melting of a finite paraffin slab in a rectangular cell. 
Problems were encountered in reproducibility of experimental 
data, due to the presence of air- bubbles in the test material 
Because of these problems, it was difficult to be certain 
that the theoretical model was actually modeling the phase- 
change phenomena. Therefore, the present study has been 
undertaken to determine whether or not the model, coupling 
of an ideal-viscous flow model with the energy equations, 
will actually predict the experimental phase change phenomena 
Also, the investigation is studying the phase change problem 
in another coordinate geometry, a cylindrical geometry. The 
experiments have been designed as to minimize the effects of 
air bubbles in the test cell. 

All phase-change experiments, such as ground tests made 
in high gravity fields, must take into account the effect of 
gravity-induced free convection. Either the experiments 
must be designed to eliminate convection or the convection 
must be mathematically modeled. It is important to deter- 
mine at what gravity level gravity-induced free convection 
may be neglected. This will enable designers of phase-change 
thermal control devices for spacecraft to determine whether 
or not gravity-induced free convection is an important design 
factor under low gravity conditions such as periods of thrust 
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Other effects, such as electrically-induced convec- 
tion or magnetically-induced convection, may also be import- 
ant design factors. Since experiments to study these other 
effects will be made in a high gravity field, the effect of 
gravity must be determined before effects of these other 
forces can be studied completely and accurately. 
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Literature Survey 


Many articles have appeared in the literature related 
to solid-liquid phase change problems and to gravity-induced 
free convection. This literature survey deals with only a 
small portion of the published material. The first treatment 
of liquefaction and solidification phenomena is that presented 
by Carslaw and Jaeger (3) in which they discuss the moving 
interface surface involved in the phase change problem. 
However, they developed only approximate mathematical models 
for semi-infinite bodies; no exact solutions were given for 
finite bodies. 

Two studies have previously been completed at this 
institution which concern the unidimensional phase-change 
phenomena of high-chain normal paraffins. Both studies use 
the one-dimensional interface equation given by Arpaci (4). 

The first study was completed by P. R. Pujado (5). In his 
study Pujado presented a theoretical model for the unidimen- 
sional melting of a finite paraffin slab. The theoretical 
model was developed using finite difference methods to 
approximate the solution of the partial differential equa- 
tions governing the physical system. The finite difference 
approximations were solved on an IBM-Model 360 digital com- 
puter. The model solved two-phase, unidimensional heat con- 
duction equations with a moving interface and variable thermal 
properties. Mr. Pujado stated that the theoretical model 
neglected free convection in the- liquid phase portion of the 
system and concluded that the errors in his results were 
probably due to the existence of free convection in the test 
cell. Mr. Pujado compared the results of his study to an 
investigation conducted earlier by Northrop Corporation (1) 
and found the results agreed very closely with the earlier 
work. The second study was performed by Ukanwa, Stermole, 
and Golden (6). The investigation concerned the solidifica- 
tion of a finite amount of liquid paraffin. A unidimensional 
model was established for the solidification of the liquid 
paraffin, based on the numerical solution by computer of the 
two-phase heat-conduction equations with moving interface 
and movable boundary conditions. Constant properties were 
assumed for each phase. The model neglected gravity- 
induced convection, supercooling, and nucleation effects. 

A comparison was made between theoretical and experimental 
temperature profiles. Good agreement was obtained between 
theory and data, although the numerical results indicate a 
faster rate of solidification than that observed experiment- 
ally. 


A study by Shah (7) investigated solid-liquid phase 
change using microphotographic equipment and temperature 
response data as analysis tools. A two-dimensional mathe- 
matical model was developed for temperature response of the 
test cell and the average interfacial velocity during the 
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solidification process. Because no interface equations have 
been developed an approximate solution was used to calculate 
the phase change energy change. A presentation of various 
phase change calculations is given by Dusinberre (8). Com- 
parison of the results obtained from the theoretical model 
and experimental results for the temperature response 
yielded reasonably good agreement. The temperature of the 
interface predicted by the theoretical model was always 
slightly higher than the experimental data. 

Grodzka and Fan (9) listed several areas of study when 
attempting to solve the problem of free convection in phase 
change thermal control equipment. They stated that free 
convection might be induced through the following forces: 
gravity, surface tension, electricity, or magnetism. 

Some of the texts which are good theoretical references 
for convection are Carslaw and Jaeger (3), Schlichting (10), 
and Longwell (11). Bird, Stewart, and Lightfoot (12) was 
used as the reference for free convection between infinite 
parallel plates. Vallentine (13) was used as the basic ideal 
flow reference for the development of the ideal-viscous flow 
model. The majority of work on free convection effects in 
liquids and gases has been done for infinite plate systems. 
Models for this type of system have been developed by Bodoia 
and Osterle (14), Dropkin and Globe (15), Dropkin and 
Somerscales (16), Gebhart (17), Koh and Price (18), and 
Samuels and Churchill (19). 

Various papers have also been published which deal 
with the melting of finite slabs. Chi-Tien and Yin-Chao 
Yen (20) developed approximate theoretical solutions for 
temperature distributions and melting rate when the mode of 
heat transfer was natural convection caused by buoyancy 
forces. They gave numerical solutions for various ice-water 
systems. Goodman and Shea (21) used a series solution to 
solve the problem of unidimensional melting of a finite slab. 

Wilkes and Churchill (22) made a study of temperatures 
in a closed rectangular system to determine the theoretical 
effects of gravity induced convection. The theoretical 
model was developed from the basic equations of motion, 
energy, and continuity; the assumption of a two-dimensional 
flow pattern precluded the study of turbulent flow. The 
system of equations was solved by an implicit alternating- 
direction technique . developed by Peaceman and Rachford (23). 
Instabilities in the numerical solution were noticed above 
certain Grashoff numbers. Fromm (24) developed another 
finite difference formulation for flow induced by a moving 
wall . 


Papers have also been published which discussed other 
causes of free convection besides gravity. Emery (25) has 
studied magnetically induced convection. Pearson (26) and 
Nield (27) have studied the effects of interfacial tension on 
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convection. They concluded that for a gas-liquid interface 
surface tension was the controlling factor for thin liquid 
layers. At some critical liquid layer thickness gravity 
induced convection becomes the controlling factor. This may 
also be true for a solid-liquid interface, but no study of 
this phenomena has been reported in the literature. 

An earlier study (2) performed at Colorado School of 
Mines dealt with the problem of gravity-induced free convec- 
tion in the melting of a finite paraffin slab with cartesian 
geometry. An ideal-viscous flow model was assumed to model 
the flow pattern; this flow model. was coupled with the energy 
equations, in finite difference form, to give a theoretical 
solution. Problems were encountered in reproducibility of 
experimental data, du-e to the presence of air bubbles in the 
test material. Because of these problems, it was difficult 
to be certain that the theoretical solution was modeling 
the phase-change phenomena accurately. 

A. 0. Ukanwa (28) performed a study to determine the 
effect of gravity-induced free convection upon the solidifica- 
tion of a finite paraffin slab. The mathematical solution 
coupled an assumed flow pattern, modified by gravity level, 
and the equations of motion into a finite-difference approxi- 
mation. Close agreement was obtained between theory and 
experimental data. A pseudo-heat capacity was used to cal- 
culate the change of phase. In personal communication with 
Dr. Ukanwa, he has indicated that the magnitude of the energy 
change involved in the phase change was modified in the com- 
puter solution to correct for impurities in the test 
material and to correct for a solid-solid phase transition. 

Lanz (29,30) and Von Rosenburg (31) have presented the 
results of studies which show that a numerical solution of a 
partial differential equation which includes bulk flow terms 
will include inherent numerical dispersion effects. These 
effects are caused by neglecting the second order partial 
differential in the finite difference approximation for the 
partial differential of time. 
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Theory 


In this study a cylindrical-coordinate finite-difference 
model is developed to predict the transient temperature 
response of a phase change material to a step change in one 
boundary temperature. The model presents one method of cal- 
culating the effect of gravity-induced free convection upon 
the phase change process. 

The test material used in the phase change investiga- 
tion was n-octadecane . The literature values for the physical 
properties are taken from an earlier study ( 2 ). 

Density 

Solid phase = (-0 . OOO 8336 )T + 1.0918, gm/cc 

Liquid phase = (-0 . 0012505 )T + 1.1316, gm/cc 

Heat Capacity 

Solid phase = 2.164, Watt-sec/gm °K 

Liquid phase = (0.008213)T - 0.14237, Watt-sec/gm °K 

Thermal Conductivity 

Solid phase = (-0.50054 x 10~^)T + .002914 watt/cm °K 

Liquid phase =(-0.50054 x 10 - ^)T + ,002914 watt/cm °K 

Melting Point = 300.60 °K 

Liquefaction enthalpy = 243.893 watt-sec/gm 

A diagram of the nodal system and boundary conditions is given 
in figure 1 . 

Ideal-Viscous Flow Model 

An earlier study (2) has shown that it is impractical 
to solve the analytical equations of motion governing the 
liquid phase when the driving force is gravity. Therefore, 
an ideal-viscous flow model is developed to approximate the 
flow pattern which exists in the liquid phase. The velocity 
profile used in the model is an approximate profile obtained 
by combining an ideal flow system for flow in a cul-de-sac 
region ( 13 ) with a viscous flow solution for flow between 
infinite parallel plates. A maximum velocity is imposed on 
the ideal-viscous flow pattern, using either a driving 
velocity calculated from the buoyancy force term or the 
maximum velocity calculated from the liquid phase energy 
equation stability criteria. 



Heating Chamber 



Figure 1. Model Network Diagram 
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The velocity profile for unidimensional flow between 
infinite vertical plates is developed by Bird, Stewart, and 
Lightfoot (12). The velocity profile is 


v 


P 


p$gb 2 AT 

12p 


(n 3 -n) 


where n = y/b . 


( 1 ) 


Since a maximum velocity has been determined the velocity 
profile should be given as a function of the maximum velocity. 
The equation for velocity now becomes 


v 


P 


v 

max 



(2) 


where n = - — (from centerline) 

° /3 

The ideal flow velocity profiles are developed by the 
use of complex variable transformations; the flow pattern 
under consideration is shown in the following diagram. 


B A 



=C1 

=C2 

=C3 
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The basic assumptions made in the conformal transforma- 
tion process are that (1) the flow pattern being studied is 
irrotational flow of a perfect fluid and that (2) the com- 
plicated flow pattern can be transformed by use of complex 
variables into parallel uniform flow. 

The flow pattern shown above is assumed to be a complex 
z-plane flow within straight-wall boundaries. Since we are 
investigating ideal flow within a simple polygon, a Schwarz- 
Christoffen transformation (13) may be used to obtain a 
parallel uniform flow pattern. 

If the polygon is in the z-plane and the new plane is 
the t-plane, then 

S§ - A 1 (a~t )”' ,r/lT (b-t )*' e/lr Cc-fc ) - ^ /lT (3) 

where A* = a complex constant 

a,b,c = real constants in ascending order of magnitude 
Y,e,£ = external defection angles of the polygon 

for the flow pattern under study 

!_=_£=_£=_! 

7T IT TT ” 2 

The boundary conditions for the transformation are 


g A, 

t 

= — oo 

g B, 

t 

= -1 

g c, 

t 

= +1 

g D, 

t 

= +oo 

Therefore , 




z = A' / — — — : + B’ (4) 

/(1-t) (-1-t) 

or z = A’ cosh -1 (t) + B’ (5) 

Applying the boundary conditions 
g C, t = 1, z = 0 
0 = cosh -1 ( 1 ) + B 1 
. * . B’ =0 
g B, t = -1, z = il 
il = A' cosh -1 (-l) 
il = A'Itt 
A’ = 1/ TT 

Therefore , 

z = -~p cosh -1 (t) 


( 6 ) 



10 


According to definition, as given by Vallentine (13), 
uniform ideal flow should have a source at and a sink at 
+°°. However, the above t-plane flow pattern gives a source 
at +°° and a sink at -<=°. Therefore, w = -t, where w is a n ew 
complex plane, and 

w =-cosh(^) (7) 

where z = x + iy 

( 8 ) • 

(9) 

( 10 ) • 

ip = -sinh sin (11) 

The stream function is defined by the following equations 

u = iP y (12) 

v = -1> x (13) 

Substituting into equation (11) and differentiating, we obtain 
u = - y sinh (~) cos(p-) (1*0 

and v = y cosh(pO sin(p-) (15) 

The liquid phase is split into three regions, as shown below. 


then 

. . ttx , iiry . 
w = -cosh(Y~ + — ~~ ) 

or w = -cosh(p-)cos (— )-i sinh(=p)sin(jpO 

By definition of complex flow 
4> = -cosh (^) cos (p^) 


REGION 

REGION 

REGION 

I 

II 

III 



► 




z 




'r 
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Region I flow is governed by equations (11), (1*0, and 
(15) with appropriate boundary equations. Region II is 
governed by equation (2) with a given maximum velocity. 

Region III is governed by another set of ideal flow equations; 
but assuming symmetrical flow it is not necessary to develop 
the equations for this region. 

The ideal flow regions are coupled to the viscous flow 
region by assuming the velocities in the viscous flow region 
are the boundary values for the ideal flow region. All y 
velocities are zero at this boundary. By use of equation (11) 
values of the stream function may be calculated at the ideal- 
viscous flow boundary. Since there can be no flow across a 
line of constant stream function, velocities in the ideal flow 
region can be related to boundary velocities at calculated 
values of the stream function at the boundary. By this 
method a pseudo-viscous flow pattern can be imposed upon the 
ideal flow regions. 


The actual liquid phase in the test cell does not have a 
constant depth. The velocities calculated by the ideal- 
viscous flow model are imposed on the liquid phase at any 
point by assuming the depth of the liquid at that point to be 
the depth of the cell in the ideal-viscous flow model. 

The flow pattern calculated is only an approximation, 
but with the small magnitude of allowable velocities calcu- 
lated by stability criteria the velocities calculated should 
give fairly accurate flow patterns in the liquid phase. 


Finite Difference Formulation of the Energy Equations 


The basic energy equations (12) governing the solid and 
liquid phases are given below. 


Solid phase 


T, = a (T + — T + T ) 
t s zz r r rr 

Liquid phase 

T. + uT + vT = a, (T + - T + T ) 
t z r 1 zz r r rr 

Boundary conditions 

@ z = 0, q| z = 0 


(16) 

(17) 


z = L, q 


= 0 


% r = R, T = T 
§ r = d , T = T 


r = r 2 > q lr 


f 

= 0 


Because of the small node size needed in the finite 
difference solution, an implicit method of solution is 
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developed. The implicit formulation eliminates the stability 
criteria with respect to time step. An implicit alternating- 
direction technique has been developed by Peaceman and 
Rachford (23) which is inherently stable with respect to time 
and spatial increments. 


The method involves the use of two successive time steps, 
each of duration At. Over the first time step, the deriva- 
tives in the R direction are solved implicitly, and the 
derivatives in the z-direction are solved explicitly. The 
procedure is reversed for the second time step. Let T* 
denote temperatures computed at the end of the first time 
step, and T° denote values at the end of the second time step. 
The r direction has an j subscript; the z direction has an i 
subscript. The following definitions are needed: 

i = 1,N; where 1 and N are boundaries, Nj_ = l s ^ liquid 
node in z direction for any given j . 

j = 1,M; where = 1 st solid node in r direction for 
any given i, and 1 and M are boundaries. 

The solid phase finite difference approximations of 
equation ( 1 6 ) are given below. 

First Time Step (Implicit in r - explicit in z) 

§ an interior node 


A 1 T*(i,j-l) + B 1 T*(i,j) + C^dJ+l) = D 1 


where 


a At , 1 

A, = -f — (•=— + i— ) 
1 Ar Ar 2r 


B 1 = 1.0 + 


2a At 

* 


a At , -i 

n — _ s ( i — i — ) 

U 1 Ar l Ar 2r' 


(18) 

(19) 

( 20 ) 

( 21 ) 


a At 

D, = T(i,j) + (T( i+1 , j ) + T( i-1 , j ) -2T( i , j ) ) (22) 

1 ( Az ; 


@ j = M. 

J 


B 1 T*(i,j) + C 1 T*(i,j+l) = D 2 

(23) 

where 


d 2 = D i - A 1 T(i,j-l) 

(2H) 

@ j = M-l 


A 1 T*(i,j-l) + Bi T* (i , j ) = D 3 

(25) 


where 
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D3 = D i - C^U.M) 


(26) 


If only one solid node is present for any given i, then 

T*(i,M-l) = D^, j = M-l (27) 

where 

= (D 3 - A 1 T(i,M-2))/B 1 (28) 

If only two solid nodes are present for any given i, then 


where 


T*(i,j) = F ± /G 1 , Dor 

j = M-l 


(29) 

T* ( i , j ) = D 5 /B 1 -C 1 T*(i 

,M-1)/B 1 , 

for j = M-2 

(30) 

Dg = D i - A 1 T(i,j-l) , 

for j = 

M-2 

(3D 

D 6 = D 1 - C 1 T(i,j+l) , 

for j = 

M-l 

(32) 

F 1 = °5 /B l - °6 /A l 



(33) 

G 1 = C l /B l " V A 1 



(34) 


Second Time Step (implicit in z - explicit in r) 
§ an interior node 


A ? T° ( i-1 , j ) + B ? T° ( i , j ) + C 7 T°(i+l,j) = D ? 


where 


A 7 


a At 
s 

(Az)‘ 


b 7 = 1.0 + 


2a At 

* 


C 7 


a At 
s 

( Az ) ‘ 


a At 

D = -2 (T* (i , j+1 )+T* ( i , j-1 )-2T* ( i j j ) ) 

1 ( Ar ; 

a At 

+ 2FAF (T* ( i ,j+l )-T* (i , j-1 ) ) +T*(i,j) 

% i = 2 

B ? T° ( i j j ) . + C 7 T°(i+l,j) = Dg 


(35) 

(36) 

(37) 

(38) 


(39) 


( 40 ) 


Dg = D ? - A ? T* ( i-1 , j ) 


where 


(41) 
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@ i = N. - 1 

A 7 T°(i-l,j) + B 7 T°(i,j) = D 9 (42) 

where 

D 9 = D ? - C 7 T*(i+l,j) (43) 

If only one solid node is present for any given i, then 

T° ( 2 , j ) = D 1() , i = 2 (44) 

where 

D 1q = (D ? - A 7 T*(l J j))/B 7 (43) 

If only two solid nodes are present for any given j , then 

T° (i, j ) = F ? /G 7 , for i = 3 (46) 

T°(i,j) = D 13 /B - C ? T° ( 3 , j )/B ? , i = 2 (47) 

where 


11 

= D ? - A y T* ( i-1 ,j ) , 

for i = 2 

(48) 

12 

= D ? - C y T* ( i+1 , j ) , 

for i = 3 

(49) 

7 

D ll /B 7 " D 12 /A 7 


(50) 

7 

c ? /b 7 - b ? /a 7 


(51) 


In the computer program, all coefficients that have r in the 
denominator must be calculated as a function of radial posi- 
tion. 


The liquid phase finite difference approximations of 
equation (17) are given below. 

First Time Step (implicit in r - explicit in z) 

§ an interior node 


A 13 T*(i,j-l) + B 13 T*(i,j) + C 13 T*(i,j+l) = D 


13 


where 


A = ^ f. v (i>j ) 

13 Ar v 2 


a . a 

— + — ) 
Ar 2r ' 


b 13 = 1.0 + 


2a 1 At 

(Ar) 2 


r - At ( zLi-iAl 

13 Ar 1 2. Ar 2r ; 

D 13 = A ~ 2Az ~ " J ' ~ CTU-i.j) - T(i+1, j ) ) + T ( i , j ) 
Ata 

+ ~(T(i+l,j) + T( i-1 , j ) - 2T( i , j ) ) 

(Az) 2 


(52) 

(53) 

(54) 

(55) 

(56) 
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@ j = 2 


where 


B 13 T*(i 5 j) + C 13 T*(i,j+l) = D l4 


D 14 = D 13 " A 13 T(i ’ 1} 


J - M - 1 


A 13 T«(i,j-l) + B 13 T*(i,j) = D 


where 


D 15 ' D 13 - A 13 T(1 >V 

If only one liquid node id present for any given i, then 
T*(i,2) = D l6 , j = 2 


where 


= (D, q - C n ,T(i,3))/B. 


(57) 

(58) 

(59) 

(60) 
(61) 
( 62 ) 


16 '~15 

If only two liquid nodes are present for any given i, then 


where 


T*(i 

,j) = F 13 /G 13 , for j = 3 

(63) 

T* ( i 

,j) = D 17 /B 13 - C 13 T»(i,3)/B 13 , j = 2 

(64) 

°17 

= D 13 - A 13 T(i,j-l), for j =2 

(65) 

co 

i — i 

Q 

= D 13 - C 13 T(i J j+l) i for j = 3 

(66) 

F 13 

D 17 /B 13 " D l8 /A 13 

(67) 

G 13 

G 13 /B 13 " B l 3 /A 13 

(68) 


Second Time Step (implicit in z - explicit in r) 
§ an interior node 


where 


A 19 T°(i-l,j) + B 19 T°(-i,j) + C 19 T°(i+l,j) = D 19 


A = At , u(i,J ) !lv 
19 hz 2 ~ Az J 


b 19 = 1.0 - 


2a^At 
(Az) 2 


(69) 

(70) 

(71) 


p -At , u ( i , j ) 
°19 Az 1 2 


- — ) 

Az' 


( 72 ) 
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D 19 = Y ' ( 2 ~ Ar )At - T*(i,j+D) + T*(l,j ) 

a. At 

+ — p (T*(l,j+1). + T* ( i , j -1 ) - 2T*(i,J)) 

(Ar)^ 


a. At 



a, at 

+ 2FAF (T* ( i , j +1 ) - T*(i,j-1)) 

(73) 

§ i = 

N i 



B ig T°(i,j) + C ig T°(i+l,j) = D 20 

(74) 

where 


D 20 = D 19 " A 19™ i ' ljJ) ' 

(75) 

@ i = 

N-l 



A ig T°(i-l,j) + B ig T°(i,j) = D 21 

(76) 

where 


°21 D 19 c ig T5t ^ 1+1 »J ^ 

(77) 

If only one liquid node is present for any given j , 

then 


T°(i,j) = D 22 , i = N-l 

(78) 


where 

°22 = (D 21 ‘ C 19 T * (N ’J)) /B i 9 (79) 

If only two liquid nodes are present for any given j , then 

T° (i , j ) = F ig /G 195 for i = N-2 (80) 

T°(i,J) = D 23 /B 19 - C 19 T°(N-2,j)/B 19 , i = N-l (8l) 

where 


23 

= D ig - A ig T*(i-l,j), i = N-2 

(82) 

24 

= D ig - C ig T*(i+l,j), i = N-l 

(83) 

19 

D 23 /B 19 " D 24 /A 19 

(84) 

19 

C 19 /B 19 " B 19 /A 19 

(85) 


Consider the solid phase energy equation. In particular, 
let equation 18 be applied to each point j = 3, 4 , . .., M-2 
in the i th column. Equations (23) and (25) hold for j = 2 and 
j = M-l, respectively. A system of simultaneous equations 
results, but there are a maximum of three unknowns in each 
equation. A non-iterative method for the solution of this 
system of equations, known as the tridiagonal matrix regres- 
sion technique, is available. These equations are of the form 
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b 2 S 2 + C 2 S 3 


d 2 

a 3 S 2 + b 3 S 3 + C 3 S 4 


- d 3 

• • 

a M-2 S M-3 + b M-2 S M-2 

• 

+ C M-2 S M-1 

d M-2 

a M-l S M-2 

+ b M-l S M-l 

d M-l 


( 86 ) 


where s. = unknown temperature at node i,j. The values of a., 
bj , and^cj are determined from equations (19), (20), and (21^. 
The values of dq , d j , and djv[_q are determined from equations 
(2*1), (22), and (26; respectively. The matrix of the coef- 
ficients of temperature is a trigonal matrix. The solution of 
equation (36) takes advantage of the tridiagonal properties of 
the coefficient matrix. The values of sj satisfying equation 
(86) are given by 

S M-1 = g M-l 

S 1 = g l " f l s l+l 5 for 1 = M ~ 2 ’ M_3 > *> 3,2 (88) 


where the g's and f's are determined by the recursion formulae 
w 2 = b 2 (89) 

W 1 = b l " a l f £-l’ for 1 = 3,4, . . .M-2,M-1 (90) 

f l = c i /w i’ for 1 = 2,3, . . . ,M-2,M-1 (91) 

g 2 = d 2 /w 2 ( 92 ) 

S 1 = (d l" a l g l-l )/w l 5 for 1 = 2 ,3,. ..,M-2,M-1 (93) 


This regression technique is used to solve both solid 
and liquid phase energy equations. During the first time step 
each row, j = 2,...,M-1, is solved by the above technique, 
going explicitly from i = 2,...,N-1. During the second time 
step the procedure is reversed, using the appropriate energy 
equations . 


This solution of the energy equations, coupled with the 
phase change calculations, gives a theoretical solution of 
the liquefaction of the test material under the influence of 
gravity-induced free convection. 

Phase Change Calculations 


A summary of various phase change calculation techniques 
is given by Dusinberre (8). A variation of the method of 
"excess degrees" was used as a calculation procedure in this 
study. Since the n-octadecane used in the study was practical 
grade, and not the pure material, the assumption was made that 
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the phase-change material changed over a 1.76°K degree tem- 
perature range. The phase-change temperature range was made 
symmetrically around the literature value for phase change 
temperature. For n-octadecane the heat capacity is the same 
above and below the phase-change temperature; in this case 
when the latent heat is divided by the heat capacity we get 
a term with the dimensions of temperature, called "excess 
degrees." It is the temperature rise which would occur if 
the amount of heat, equal to the. latent heat, were added and 
no phase change took place. The procedure for calculating 
the phase change is given below. After each iteration a test 


is run on solid phase temperatures. 

Given: 

T s (i,j).R.T f (94) 

If, T s (i,j) < T^, the node is still solid (95) 

If, T s (i, j) >. T~, the node is changing phase (96) 

If equation (96) is applicable for the node being 
investigated, then the following procedure is followed. 

Given: 

Tg ( i , J ) * R 4 AH f /C p (97) 

If, T (i,j) < AH f /C , the node has not changed phase (98) 
and the S temperature is given by 

T s (i,J) = T fo + T e (i,j)*C p *1.76/AH f (99) 

If, T e (i,J) > AH f /C p (100) 

the node has changed phase 'and the liquid phase temperature 
is given by 

T^iJ) = T fo + 1.76 + (T e (i,j)*C - AH f )/C p (101) 

Numerical Dispersion Effects 


Von Rosenburg (31) and Lantz (30) state that finite 
difference solutions of partial differential equations which 
have bulk flow velocity terms include numerical dispersion 
effects caused by ignoring the second order partial differ- 
ential with respect to time in the Taylor’s series expansion 
used to derive the finite difference approximation for the 
partial of the variable with respect to time. Both explicit 
and implicit solutions exhibit these numerical dispersion 
effects. 

When the second partial with respect to time is left 
in the Taylor’s series expansion and a partial differential 
equation is derived from the finite difference equation, 
equation (17) becomes 
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. , a.. 

T, + - T, , + uT + vT = — T + a, T + a.,T (102) 

t 2 tt z r r r - 1 rr 1 zz 


The procedure for calculating the numerical dispersion terms 
is 

(1) rearrange equation (17) and solve for 

(2) substitute T t , from step 1, into equation (102) 
for term — 


(3) rearrange the resulting equation. 
When this is done, equation (102) becomes 


A-t- on At a. 

T, + uT + (v + — )T = — (1.0 + — ~) T 

t z 2r ■ r r „ 2 r 

2r 

O 

Atonv .. 2 Aton^ 

+ (a + L_ - — ) T 

'‘l 2r 2 2r ;i rr 


+ a l(10 - ^)T zz 


Atva 


+ (■ 


1 


- Atvu)T 
r rz 


(103) 


To minimize the effect of numerical dispersion the magnitude 
of the time step in the finite difference solution must be 
reduced until numerical dispersion coefficients, all coef- 
ficients which include At in the above equation, are negligible, 
and a convergent solution is obtained when the time step is 
varied. 


Stability Criteria 

Although the implicit solution developed in'/this study 
has eliminated the time step stability criteria, stability 
criteria still exist which limit the magnitudes of allowable 
velocities (30) in the program. The stability criteria for 
velocities from the finite difference formulation of equation 
(17) are given, below. 


2a^/Az 

(104) 

a (r/Ar-1 ) 

, for v < 0 

(105) 

a (r/Ar+1) 

, for v > 0 

r» ? 

( 106 ) 
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Equipment and Procedure 


In this section a short discussion of the equipment 
and procedure used in the experimental investigation portion 
of this study is given. The test cell, figures 2 and 3, 
consisted of an annular test chamber, a tubular heating 
chamber, and an expansion chamber. The heating chamber was 
aluminum pipe with an outside diameter of 1.906 cm. The out- 
side diameter of the test chamber was 5.08 cm; the height of 
test chamber was 5-08 cm. The test chamber was drilled in a 
10.16 cm by 10.16 cm by 7.62 cm block of plexiglas. The 
expansion chamber, 7-62 cm by 7.62 cm by 2.54 cm, was machined 
from a 10.16 cm by 10.16 cm by 3-76 cm block. Sixteen iron- 
constantan thermocouples, made from 24-gauge wire, were placed 
in the test chamber at various positions. These positions are 
given below. The z-position is the distance from the bottom 
of the test chamber; the r-position is the distance from the 
inside wall. 


Thermocouple 

z-position 

r-position 

No . 

(cm) 

(cm) 

1 

1.016 

1.270 

2 

1.016 

0.9525 

3 

1.016 

0.635 

4 

1.016 

0.3175 

5 

2.032 

1.27 

6 

2.032 

0.9525 

7 

2.032 

0.635 

8 

2.032 

0.3175 

9 

3.048 

1.27 

10 

3.048 

0.9525 

11 

3.048 

0.635 

12 

3.048 

0.3175 

. 13 

4.064 

1.27 

14 

4 . 064 

0.9525 

15 

4.064 

0.635 

16 

4.064 

0.3175 


Thermocouples were also used to measure the heating tank tem- 
perature and the heating chamber wall temperature. The 
expansion chamber and test chamber were connected using 
4-0.3175 cm diameter bolts. A cork gasket was placed in 
the test chamber to ensure that no leaks developed. Epoxy 
was used to seal all other openings, such as the thermocouple 
ports and the heating chamber port. 

The heating system consisted of the following equip- 
ment: a constant temperature bath, a constant temperature 

controller, a centrifugal pump, a flow-meter, and lines and 
valves. A diagram of the heating system is given in figure 4. 
The flow lines consisted of 1.27-cm inside-diameter copper 
tubing; the valves were 1.27 cm Prier globe valves. The flow- 





5.08 cm dia 


4-0.3175 cm 
dia bolts 


Figure 3. Test Cell - Top View 





Figure 4. Flow System 
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meter was a Fisher and Porter Co. precision bore flowrator 
tube No. FP-1/2-27-G-10/83 with a 100% rated capacity of 
2.3^67 liters/min for a liquid with a specific gravity of 

1 . 0 . 


The constant temperature bath was a 5-liter Pyrex 
tank. The constant temperature bath was a Hanke Company, ' ' 

Model E51, Constant Temperature Circulator. The circulator 
is designed for thermosetting open baths in the temperature 
range 243 to 423°K with a control accuracy of ±0.02°K. The 
circulator has a continuously variable heater output selec- 
tions between 0 and 1000 watts. The circulating pump used to 
circulate the heating fluid to the test cell was a Chemical 
Rubber Company "No Seal" centrifugal pump. Model AB1P005N. 

The pump operated on 115-volt, 60 cycle/min, alternating 
current. The pump's rated capacity was 2.65 liters/min at a 
head of 30.5 cm to 1.58 liters/min at 274.5 cm under normal 
atmospheric conditions. 

Temperatures were recorded using a Bristol Dynamaster 
Multipoint convertible recorder. Model 570, operating on 120 
volt, 60 cycle/min, alternating current. The print speed 
was 2 seconds per thermocouple. The accuracy of the recorder 
was ±0.417 °K. 

Experimental runs were made using the following procedure 

1. For a melting run the tank temperature was set 
approximately 1°K higher than the desired hot wall tempera- 
ture. The pump was turned on and heating fluid, water, was 
allowed to flow through the bypass lines; this brought most 
of the water in the flow system to tank temperature. When 
the constant temperature bath reached a constant temperature, 
and the test cell thermocouples recorded temperatures that 
were within 0.278°K the experimental run was started. For a 
solidification run the test material was heated 5*6°K above 
the melt point. The constant temperature bath was filled 
with ice, and the flow cooled to 273.l6°K. When the test 
material cooled to 2.78°K above the melt point the experi- 
mental run was started. 

2. At the start of the experimental run the test cell 
was first leveled; then the run was started by switching from 
bypass flow to heating chamber flow. The starting point was 
marked on the chart paper of the temperature recorder. 

3. The duration of the solidification runs was 
approximately forty minutes; the duration of the liquefac- 
tion runs varied from 50 minutes to 120 minutes. 
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Discussion of Results 


In this section the results of two experimental solidi- 
fication runs and ten experimental liquefaction runs are 
compared to theoretical temperature profiles obtained using 
the numerical computer solution. The experimental solidifi- 
cation runs were made to determine a latent heat of fusion 
under pure conduction conditions. All liquefaction experi- 
mental runs show the effect of gravity-induced free convection. 
Good agreement was obtained between theory and data for the 
liquefaction runs; agreement was best when the liquid-phase 
temperature gradient was small. 

Figure 5 shows typical results from the combined experi- 
mental and theoretical determination of a latent heat of 
fusion to be used in the liquefaction study. During personal 
communication with A. 0. Ukanwa and S. P. White ( 32 ) it was deter 
mined that the use of literature values of latent heat in the 
numerical solution of solid-liquid phase change problems 
being studied was not modeling the physical situation being 
observed experimentally. In all three investigations, two 
liquefaction studies and one solidification study, the test 
materials were all high chain normal paraffins of. practical 
chemical grade. There are three possible justifications for 
varying the literature value of latent heat; (1) Since only 
practical chemical grade test material was used in the study, 
the physical properties may be different from literature 
values. (2) Impurities may exist in the test material, due 
to chemical reaction with the aluminum walls, leaching of 
solvents from the plexiglas walls, or other forms of contamina- 
tion. (3) In an earlier study (2) and in the present study, 
a large number of air bubbles have been observed in the test 
material during the experimental runs. These air bubbles 
have an effect on the physical properties of a given volume 
of the test material. 

In the computer solution the temperature response is 
insensitive to thermal diffusivity, making the latent heat 
of fusion the governing factor. The results of the compari- 
son of theory and data shown in figure 5 indicate that a 
latent heat of 75 percent of literature value gives the best 
theoretical approximation of experimental data. Ukanwa (28) 
overcame the problem by using a pseudo-heat capacity to model 
the phase change; he stated that the method, as used in the 
computer solution, also changed the magnitude of the latent 
heat. However,, this method of phase change calculation is 
not applicable to a liquefaction study. White, using n- 
octadecane as a test material, has also obtained good agree- 
ment between theory and data using a latent heat of 75 per- 
cent of literature value. The only deviations that appear 
between theory and data, figures 5-a, 5-b, and 5-c, occur at 
the phase change temperature; this indicates that a larger 
temperature range for phase change may be applicable. Until 
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the phase change takes place there is little apparent dif- 
ference in the theoretical profiles for the different latent 
heats, see figure 5-d. Only one thermocouple is presented 
for each radial position; since all thermocouples for both 
experimental runs were within 0.7°K at any radial position 
the solidification was unidimensional, and only one thermo- 
couple was needed in the presentation. The experimental • 
data, runs C-19-2 and C-23-2, are given in Appendix A. As in 
the equipment and procedure section the z-position of a 
thermocouple is the distance from the bottom of the test 
cell and the r-position is the distance from the inside wall. 

Figure 6 presents a comparison of experimental tempera- 
ture data and theoretical temperature predictions for a 
liquefaction run with a hot wall temperature of 313.55°K. 

The agreement between experimental and theoretical results 
is very good for all sixteen thermocouples. Good reproduci- 
bility of experimental data is shown by the two experimental 
runs presented. The effect of convection is very important. 
For example, figures 6-a, 6-b, 6-c , and 6-d are all at 
R = 0.3175 cm and various z-positions. If the mode of heat 
transfer were pure conduction then all four thermocouples 
would show that the interface was flat. But thermocouple 
16 shows that the node at z = 4.064 cm melts at 1980 seconds; 
thermocouple 12 shows that the node at z = 3-048 cm melts at 
3000 seconds; thermocouple 8 shows that the node at z = 2.032 
cm melts at 4800 seconds; and thermocouple 4 shows that the 
node at z = 1.016 cm melts at 5400 seconds. The melt pattern 
is definitely affected- by convection. The thermocouples at 
other R-positions show the same effect, but with a time lag 
caused by their larger distances from the hot wall. At any 
given z-position the final liquid phase temperatures for all 
R-positions are approximately the same, see figures 6-d, 6-h, 
6-1, and 6-p. 

Figure 7 shows the effect of numerical dispersion upon 
the theoretical temperature profiles for a hot wall tempera- 
ture of 313-55°K. For the duration of the runs, the largest 
effect is at z = 4.064 cm, figures 7-d, 7-h, 7-1, 7-p. The 
results show that the solution obtained is not yet in the 
convergent region. Due to computer limitations it was not 
possible to obtain a convergent time step. With a larger and 
faster computer solution it should be possible to reduce the 
time step far enough to eliminate the effect of numerical 
dispersion. At z = 2.032 cm and z = 3-048 cm, for example 
see figures 7-b and 7-c, there is very little effect due to 
numerical dispersion, due to the fact that these z-positions 
are in the parallel flow region and not influenced by the 
velocity to a large extent. At z = 1.016 cm the solution is 
not affected by numerical dispersion, due to the fact that 
temperature gradients in this region of the test cell are . 
small. In all liquefaction computer runs the time step used 
was the smallest one that would allow modeling of an entire 
experimental run; core limitations were placed on the com- 
puter solutions that made it necessary to complete each run 
within a seven hour time limit. 
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Figure 8 presents a comparison of theoretical temperature 
profiles and experimental temperature profiles for a hot wall 
temperature of 319. 11°K. There is good agreement between 
theory and data for all thermocouples; but larger deviations 
occur in the final liquid phase temperatures than occurred at 
a hot wall temperature of 313.55°K. The largest deviations 
in final liquid phase temperatures are 2.25°K, see figures 
8-d, 8-h, 8-o, and 8-p; all phase change times are predicted 
within 180 seconds. The R = 0.3175 cm, z = 2.032 cm and 
z = 3.048 cm thermocouples, see figures 8-b and 8-e, show 
deviation between theory and data immediately after the 
phase change takes place; the theoretical curves show a more 
rapid temperature rise than the experimental data indicates 
should occur. This deviation is accounted for by the fact 
that the theoretical model assumes parallel flow in this por- 
tion of the liquid phase while experimentally the flow is not 
parallel in this region, because the interface is not parallel 
to the hot wall. 

Figure 9 shows a comparison of three theoretical runs 
made for a hot plate temperature of 319.11°K. As in figure 
7 the results show that the solution is not in the convergent 
region, and that numerical dispersion is still an important 
factor in the solution for the time steps used. The effect 
is of the same magnitude as that observed for figure 7. 

Figures 10, 11, and 12 present theoretical temperature 
profiles compared to experimental temperature profiles. As 
in earlier runs, the experimental reproducibility of data is 
very good. The same trends are present as in figure 8; all 
phase change times are predicted very closely; deviations are 
present in the rate of liquid phase temperature rise for 
z = 2.032 cm and z = 3-048 cm for R = 0.3175 cm. Again the 
model does predict accurately the final liquid phase tempera- 
tures at z = 4.064 cm. For a hot wall temperature of 321.33°K, 
figure 10, the largest deviation in -final liquid phase tem- 
perature is 2.8°K; for a hot wall temperature of 327-44°K, 
figure 11, the largest deviation in final liquid phase tem- 
perature is 4.4°K; for a hot wall temperature of 330.22°K, 
figure 12, the largest deviation in final liquid phase tem- 
perature is 5*6°K. 

From the results presented in this section, it can be 
concluded that gravity-induced free convection is an import- 
ant design factor in the use of passive solid-liquid phase- 
change thermal control devices. In all ground tests gravity- 
induced convection will affect results to some extent, even 
in cases where the experiments are designed to minimize the 
effects of gravity-induced convection. In cases where the 
hot wall is not horizontal one end of the test material, that 
at the highest elevation, will melt faster than predicted by 
pure conduction. With the present progress made in prediction 
of gravity-induced convection it is not possible to say at 
what gravity level convection effects may be neglected. 
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Figure 5» Comparison of Experimental Data for a Solid- 
ification Run to Theoretical Pure Conduction 
Profiles for Various Values of Latent Heat 

(a) at R = 0.3175 cm 

(b) at R = 0.635 cm 

(c) at R = 0.9525 cm 

(d) at R = 1.27 cm 


Legend : 

^used “ ^^literature 

K = 1.0 

K = 0.75 

K = 0.50 

O - Experimental Data 













Figure 6. Comparison of Experimental Data to Theoretical 
Model Temperature Profiles for a Hot Wall 
Temperature of 313»55°K 

(a) r = 0.3175 cm, z = 1.016 cm 

(b) r = 0.3175 cm, z = 2.032 cm 

(c) r = 0.3175 cm, z = 3.048 cm 

(d) r = 0.3175 cm, z = 4.064 cm 

(e) r = 0.635 cm, z = 1.016 cm 

(f) r = O.635 cm, z = 2.032 cm 

(g) r = 0.635 cm, z = 3.048 cm 

(h) r = 0.635 cm, z = 4.0 4 cm 

(i) r = 0.9525 cm, z = 1.0l6 cm 

(j) r = 0.9525 cm, z = 2.032 cm 

(k) r = 0.9525 cm, z = 3.048 cm 

(l) r = 0.9525 cm, z = 4.064 cm 

(m) r = 1.270 cm, z = 1.016 cm 

(n) r = 1.270 cm, z = 2.032 cm 

(o) r = 1.270 cm, z = 3.048 cm 

(p) r = 1.270 cm, z = 4.064 cm 

Legend: 

O = C-13-2 

Q = C-14-2 

= Theoretical Model 







Figure 6 (cont) 
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Figure 6 (cont) 
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Figure 7* Effect of Numerical Dispersion Upon Theoretical 
Temperature Profiles for a Hot Wall Temperature 
of 313.55 °K 

(a) r = 0.3175 cm, z = 1.016 cm 

(t>) r = 0.3175 cm, z = 2.032 cm 

(c) r = 0.3175 cm, z = 3.048 cm 

(d) r = 0.3175 cm, z = 4.064 cm 

(e) r = 0.635 cm, z = 1.016 cm 

(f) r = 0.635 cm, 7. = 2.032 cm 

(g) r = 0.635 cm, z = 3*048 cm 

(h) r = 0.635 cm, z = 4.064 cm 

(i) r = O .9525 cm, z = 1.016 cm 

(3) r = 0.9525 cm, z = 2.032 cm 

(k) r = 0.9525 cm, z = 3.048 cm 

(l) r = 0.9525 cm, z = 4.064 cm 

(m) r = 1.270 cm, z = 1.016 cm 

(n) r = 1.270 cm, z = 2.032 cm 

(o) r = 1.270 cm, z = 3.048 cm 

(p) r = 1.270 cm, z = 4.064 cm 

Legend: 

= 1.0 sec for time increment 

= 1.5 sec for time increment 

= 3*0 sec for time increment 



















TEMP (°K) TEMP (°K) TEMP (°K) TEMP (°K) 


































Figure 7 (cont) 
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Figure 8. Comparison of Experimental Data to Theoretical 
Model Temperature Profiles for a Hot Wall 
Temperature of 319.11 °K 

(a) r = 0.3175 cm, z = 1.0l6 cm 

(b) r = 0.3175 cm, z = 2.032 cm 

(c) r = 0.3175 cm, z = 3.048 cm 

(d) r = 0.3175 cm, z = 4.064 cm 

(e) r = 0.635 cm, z = 1.0l6 cm 

(f) r = 0.635 cm, z = 2.032 cm 

(g) r = 0.635 cm, z = 3*048 cm 

(h) r = O.635 cm, z = 4.064 cm 

(i) r = 0.9525 cm, z = 1.0l6 cm 

(j) r = 0.9525 cm, z = 2.032 cm 

(k) r = 0.9525 cm, z = 3.048 cm 

(l) r = 0.9525 cm, z = 4.064 cm 

(m) r = 1.270 cm, z = 1.016 cm 

(n) r = 1.270 cm, z = 2.032 cm 

(o) r = 1.270 cm, z = 3.032 cm 

(p) r = 1.270 cm, z = 4.064 cm 

Legend: 

O = C-21-2 


1 


g = c-22-2 

— = Theoretical Model 
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Figure 8 (cont) 
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Figure 9. Effect of Numerical Dispersion Upon 
Theoretical Temperature Profiles for 
Hot Wall Temperature of 319*11 °K 

(a) r = 0.3175 cm, z = 1.016 cm 

(b) r = 0.3175 cm, z = 2.032 cm 

(c) r = 0.3175 cm, 2 = 3*048 cm 

(d) r = 0.3175 cm, 2 = 4.064 cm 

(e) r = O.635 cm, z ~ 1.0l6 cm 

(f) r = O.635 cm, z = 2.032 cm 

(g) r = O.635 cm, z = 3.048 cm 

(h) r = 0.635 cm, z = 4.064 cm 

(i) r = 0.9525 cm, z = 1.016 cm 

( j) r = 0.9525 cm, z = 2.032 cm 

(k) r = 0.9525 cm, z = 3.048 cm 

(l) r = 0.9525 cm, 2 = 4.064 cm 

(m) r = 1.270 cm, z = 1.016 cm 

(n) r = 1.270 cm, z = 2.032 cm 

(o) r = 1.270 cm, z = 3.048 cm 

(p) r = 1.270 cm, z = 4.064 cm 

Legend: 

= 1.5 sec for time increment 

= 3*0 sec for time increment 

=6.0 sec for time increment 
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Figure 10. Comparison of Experimental Data to 

Theoretical Model Temperature Profile s 
for a Hot Wall Temperature of 321.33 °K 

(a) r = 0.3175 cm, z = 1.016 cm 

(b) r = 0.3175 cm, z = 2.032 cm 

(c) r = 0.3175 cm, z = 3*048 cm 

(d) r = 0.3175 cm, z = 4.064 cm 

(e) r = O.635 cm, z = 1.016 cm 

(f) r = O.635 cm, z = 2.032 cm 

(g) r = O.635 cm, z = 3.048 cm 

(h) r = O.635 cm, z = 4.064 cm 

(i) r = 0.9525 cm, z = 1.01 6 cm 

(j) r = 0.9525 cm, z = 2.032 cm 

(k) r = 0.9525 cm, z = 3*048 cm 

(l) r = 0.9525 cm, z = 4.064 cm 

(m) r = 1.270 cm, z = 1.016 cm 

(n) r = 1.270 cm, z = 2.032 cm 

(o) r = 1.270 cm, z = 3.048 cm 

(p) r = 1.270 cm, z = 4.064 cm 

Legend : 

O = C-ll-2 
0 = C-12-2 

— = Theoretical Model 












Figure 10 (cont) 
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Figure 10 (cont) 
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Figure 11. Comparison of Experimental Data to 

Theoretical Model Temperature Profiles 
for a Hot Wall Temperature of 327.44 °K 

(a) r = 0.3175 cm, z = 1.016 cm 

(b) r = 0.3175 cm, z = 2.032 cm 

(c) r = 0.3175 cm, z = 3*048 cm 

(d) r = 0.3175 cm, z = 4.064 cm 

(e) r = 0.6350 cm, z = 1.016 cm 

(f) r = O. 635 O cm, z = 2.032 cm 

(g) r = 0.6350 cm, z = 3*048 cm 

(h) r = 0.6350 cm, z = 4.064 cm 

(i) r = 0.9525 cm, z = 1.016 cm 

(j) r = 0.9525 cm, z = 3.048 cm 

(k) r = 0.9525 cm, z = 3.048 cm 

(l) r = 0.9525 cm, z = 4.064 cm 

(m) r = 1.270 cm, z = 1.0l6 cm 

(n) r = 1.270 cm, z = 2.032 cm 

(o) r = 1.270 cm, z = 3*048 cm 

(p) r = 1.270 cm, z = 4,064 cm 

Legend : 

O = C-ll-2 
0 = C-12-2 

- — = Theoretical Model 
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Figure 11 (cont) 
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Figure 12. Comparison of Experimental Data to 

Theoretical Model Temperature Profiles 
for a Hot Wall Temperature of 330.22 °K 

(a) r = 0.3175 cm, z = 1.016 cm 

(b) r = 0.3175 cm, z = 2.032 cm 

(c) r = 0.3175 cm, z = 3.048 cm 

(d) r = 0.3175 cm, z = 4.064 cm 

(e) r = O .635 cm, z = 1.016 cm 

(f) r = O .635 cm, z = 2.032 cm 

(g) r = 0.635 cm, z = 3*048 cm 

(h) r = 0.635 cm, z = 4.064 cm 

(i) r = 0.9525 cm, z = 1.016 cm 

(j) r = 0.9525 cm, z = 2.032 cm 

(k) r = 0.9525 cm, z = 3.048 cm 

(l) r = 0.9525 cm, z = 4.064 cm 

(m) r = 1.270 cm, z = 1.016 cm 

(n) r = 1.270 cm, z = 2.032 cm 

(o) r = 1.270 cm, z = 3.048 cm 

(p) r = 1.270 cm, z = 4.064 cm 

Legend: 

0 = C-17-2 
Q = C-18-2 

— = Theoretical Model 
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Conclusions 


The following conclusions have been made from this 

study : 

1. The determination of physical properties, especially 
the latent heat of fusion, is very critical for the 
proper modeling of the phase-change phenomena. When 
working with materials with a very low thermal dif- 
fusivity the latent heat of fusion becomes the 
governing physical property in the theoretical model- 
ing of the system. The latent heat may be affected 
by purity of material, either chemical purity as 
manufactured or impurities introduced by chemical 
reactions with the test cell, and by air bubbles in 
the material. When using high chain normal paraffins 
air bubbles are the major factor to be considered in 
the determination of an effective latent heat of 
fusion. 

2. Numerical dispersion is an important factor in the 
modeling of the solid-liquid phase-change phenomena 
when free convection is present. In this study we 
were not able to reduce the time step in the numer- 
ical solution to a small enough value to reach a 
convergent solution where numerical dispersion could 
be considered negligible. 

3. Other sources of numerical error are the constant 
maximum velocity and the limitation placed on velo- 
city by the stability criteria in the liquid phase 
energy equation. 

4. An earlier study (2) has shown that it is not pos- 
sible to determine a gravity level until an experi- 
mental investigation has been made to determine a 
critical Rayleigh number for the particular experi- 
mental system under study. However, the conclusion 
can be made that gravity-induced free convection is 
an important factor in the design of passive phase- 
change thermal control devices, especially when 
ground tests are to be made on a phase-change 
device. In a high gravity field gravity-induced 
free convection will cause the phase-change material 
to melt faster than predicted by a pure conduction 
model in certain portions of the cell. This will 
cause hot spots in the equipment whose temperature 

is being controlled, if allowance is not made for the 
increased melting rate in the design of the phase- 
change device. 
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5. The numerical study assumed that the velocity pro- 
files were symmetric, in magnitude but not sign, 

in the development of the ideal-viscous flow model. 
However, the experiment was made using a cylindrical 
geometry test cell. Nodes further away from the 
tubular heating wall contain more mass than nodes 
close to the heating wall. Therefore, the flow 
model should be modified to account for the fact 
that velocities near the hot wall should be larger 
in magnitude than nodes further from the hot wall. 

6. The experimental data is reproducible, which means 
that air bubbles in the test material did not cause 
hot or cold spots in the material. The only effect 
of air bubbles was on the latent heat of fusion. 

By making the heating wall vertical air bubbles rose 
to the top of the cell and did not affect heat 
transfer rates from the hot wall to the test mater- 
ial . 

7. The study has shown that the method of solution, an 
ideal-viscous flow model coupled with the energy 
equation, will model the phase change process when 
liquid phase temperature gradients are small. When 
the liquid phase temperature gradients become larger 
deviations appear between final liquid phase tem- 
peratures predicted theoretically and measured 
experimentally. Therefore, the theoretical solution 
should be considered an initial solution to the 
problem of gravity-induced free-convect ion effects 
in the solid-liquid phase-change phenomena. Further 
work is needed if a complete solution to the problem 
is to be made. 
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Recommendations 


The following recommendations are made as a result of 
this study: 

1. The computer solution used in this study should be 
modified to account for the fact that velocity levels 
should become smaller in magnitude further away from 
the hot wall. The modified computer solution should 
then be compared to the experimental results of this 
study to' see if the theoretical prediction is better 
than that using the present computer solution. 

2. An experimental study using tracer materials in the 
phase change system should be made to determine the 
actual shape of convection induced velocity profiles 
in liquefaction and solidification phenomena. The 
measurement equipment should be photographic or 
microphotographic equipment. 

3. Because of the problems encountered in the theoretical 
modeling of the convection phenomena, a material 
investigation study should be made to determine whether 
or not low density inert polymers, which exhibit 
solid-solid phase changes with high heats of transi- 
tion, would make feasible phase-change materials. 

This type of material could be modeled using only a 
pure conduction model. 

4. A study should be undertaken to determine effective 
thermal dif fusivities that would model the phase 
change phenomena without solving for velocity profiles. 

5. A theoretical investigation should be made to deter- 
mine the proper finite difference form of the equa- 
tion of , state for gravity to be introduced into the 
equations of motion. By this means the actual equa- 
tions governing the liquid phase could be solved. 
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Nomenclature 


Definition: 

Given that s = f(z,R,t), then the following definitions 
are true: 

3s „ _ 3s 3s _ 3 2 s _ 3 2 s 

S t 3t » z 3 z 5 S r 3r 5 S tt ..2 ’ S zz „ 2 

3 1 3 z 

3 2 s _ 3 2 s 

S rr * s rz 3r3z 

Parallel Flow Model: 

b = one-half of distance between parallel walls, cm 

p 

g = acceleration of gravity, cm /sec 
AT = temperature gradient between parallel walls, °K 
y = distance from centerline, cm 
n = y/.b 

6 = coefficient of thermal expansion, °K 1 
p = density, grams/cubic cm 
y = viscosity, gm cm -1 sec -1 
Ideal Flow Model: 

a,b,c = real constants in ascending order of magnitude 

A’,B' = complex constants 

t,w,z = complex planes 

u = x-direction velocity, cm/sec 

v = y-direction velocity, cm/sec 

x = spatial dimension in complex z plane, cm 

y = spatial dimension in complex z plane, cm 

Finite Difference Models: 

A,B,C,D,F,G = coefficients of IAD equations 

C p = heat capacity, watts sec gm - -*- °K"'*' 

AH = latent heat of fusion, watts sec gm 

2 

q = heat flux, watts/cm 
r = radial spatial dimension, cm 
T = temperature, °K 
t = time, sec 

u = velocity in z-direction, cm/sec 
v = velocity in r-direction, cm/sec 



z = longitudinal spatial dimension, cm 
a = thermal diffusivity, cm^ sec - ^ 

Subscripts : 

e, f,l,o,p,s = excess, fusion, liquid, cold wall, hot 

wall, solid, respectively 

Superscripts 

*,o = at end of l s ^ time step, at end of second time 
step, respectively 

Tridiagonal Matrix: 

a,b,c,d = coefficients of tridiagonal matrix 

f, g,w = regression coefficients 

Computer : 

AK2 = liquid phase thermal conductivity, Btu (ft sec ° 

AKP = wall thermal conductivity, Btu (ft sec °F) -1 

AKS = solid phase thermal conductivity, Btu (ft sec °F 

2 -i 

AL = liquid phase thermal diffusivity, ft sec 

2 -1 

AS = solid phase thermal diffusivity, ft sec 
CPL = liquid phase heat capacity, Btu (lb °F) _1 
CPP = wall heat capacity, Btu (lb °F) -1 
CPS = solid phase heat capacity, Btu (lb °F) _1 
DR = spatial increment, r-direction, ft 
DT = time increment, sec 
DZ = spatial increment, z-direction, ft 
HA = heat transfer coefficient, wall-atmosphere, Btu 
(ft^sec °F) _1 

I 

HF = latent heat of fusion, Btu/lb 
Ml = number of nodes in r-direction 
M2 = Ml - 1 

M3 = Ml - 2 

MR1 ,MR2 ,MR3 ,MR4 = r-positions for nodes to be printed 
MX1 ,MX2 ,MX3 ,MX4- = z-positions for nodes to be printed 
N1 = number of nodes in z-direction 
N2 = N1 - 1 

N3 = N2 - 1 

RMIN = radius of heating tube, inches 
R0L = liquid phase density, lb(ft^) -1 
R0P = wall density, lb(ft^) - ' i ' 

R0S = solid phase density, lb(ft3) - -*- 



T 

TA 

TI 

TN 

U. 

UR 

UZ 

V 

VMAX 


temperature at old time step, °F 

atmospheric temperature, °F 

initial temperature, °F 

temperature at new time step, °F 

flow model velocity for z-direction, ft sec - ' 1 ' 

liquid phase velocity for r-direction, ft sec 

liquid phase velocity for z-direction, ft sec 

flow model velocity for r-direction, ft sec - ' 1 ' 

allowable maximum velocity, ft sec" 1- 
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- APPENDIX A - Experimental Data 


Data from two solidification runs. and two liquefaction runs 
are presented in this section. For experimental data of the 
other eight experimental runs contact Dr. J.. 0. Golden, 
C.P.R.E. Department, Colorado School of Mines, Golden, Colo. 
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This program was written to solve a two-dimensional 
cylindrical-coordinate liquefaction problem of n-octadecane 
under the influence of gravity-induced free convection. 

The program was run on a DEC, Model PDP-10 computer. The 
program was written in general terms, except for the heat 
capacity term in the phase change calculation. To use the 
program the following length to effective radius must be 
observed 

N1 > 2M1 + 1 

See the nomenclature for definition of terms and units to 
be used in the program. 

Input File: 

1. For execution 4 = Input file or device 

2. First card - N1 ,N2 ,N3 ,M1 ,M2 ,M3 ,MR1 ,MR2 ,MR3 ,MR4 , 

with a (101) format) 

3. Second card - RMIN, with a (F) format . 

4. Third card - MX1,MX2 ,MX3 ,MX4 , with a (41) format 

5. Fourth card - DZ , DR , DT , TP , TA , TF , TI , HA , HF , with a 

(9F) format. TA and HA are dummy 
input variables, to be used if cal- 
culation of heat loss through walls 
is incorporated into the program 

6. Fifth card - CPS,CPL,CPP,AKS,AKL,AKP,R0S,R0L,R0P 

with a (9E) format. CPP,AKP,R0P are 
dummy input variables, to be used if 
a heat balance of the cell walls is 
incorporated into the program 

7. Sixth card - VMAX, with an (E) format 
Flag File: 

1. For execution 6 = Flag File 

2. Output - ’Input File* 

3. Input - FLAG, with an (F) format 

a. if FLAG < 10.0, stop execution 

b. if FLAG >_ 10.0, continue execution 
Output File: 

1. For execution 5 = output file, 6 = output file 

2. Time, t, in seconds at 120 second intervals 

3. Temperatures at input nodes at 120 second intervals 

4. Interface nodes at 120 second intervals 
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Sample Input Pile: 

81,80,79,26,25,24,7,11,16,21 

0.375 

17,33,49,65 

0.0020833,0.0020833,1.5,99.75,75.0,80.0,77.5,0.0,78.675 
0 . 517E+00 , 0 . 522E+00 , 0 . IE- 12 , 0 . 243E-04 , 0 . 2402E-04 , 0 . IE-12 
0.5349E+02,0. 482E+02 , 0 . IE- 12 

0.75E-03 



c 

c 


14 
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CYL I .-TO « I Ca.L C 0 0 R ' ■ I • J a T £ HO PROGRAM 

sssssssssssssss 

n I H £ V s i n \i ij (31,3 0 ) , V < 81 , 33 > » UH < , 3D # t'P <51*33) 

T ( HI , 3'* ) « T •< ?i . 3D ,3 ( ?i , 3 2 ) 

D I ”£I* S ! Or» .*.*) (51) . Ml ( .; t ) , -3 < 3 ? > 

A { 3 -5 ) , C ( 3 X ) , D ( 3D . « < 3D , Y ( 3? > 

OJC -Hi 5 , .;2( 'H , Y2(«l) , JR ( 5 ) , J>(S) 
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WRITE; (6,1-5) MYi ,MX2,HX3.MX4 
FORMAT HI) 

RE a 0 (A.ii) DZ,D-»0T,.TP,T4 iTF»TI,^A,hF 
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FORMAT <-~F> 
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F T 

, / i *. 


L ’ r 


OF.! 


•> I 

T P = ? # r i 

“ X , 


- 1 * f , 1 f r ' 


T ! 


- • 


12 

112 


F,' PEO F PF = ' , ~ » 1 RTU/U-2 ’/> 

IOFUT PHYSICAL t£‘>. , ... 

RE A 0 (0,12) C P 0 , ?PL » 1' P P , a \ S , a K L , a •< F 

WRITE (f.u 2) 

• FORMAT (VP 

r *; L •■! IT f -I y . * 


0. ; -0L 


1 

2 

3 


Ur" 

F,' •-Ti.:/(LF#DEG 
a<i. = •, 
LS/FT**3 ' ' 

REAP; ( 4 , 1 A ) V • ; A X 
FORMAT (F) 

CALL VI VV ( V'AX ) 

I : ' v - 1 T I ALIZE VARIABLES 
j R ( 1 ) sroi 
JR (2 )=?•!* 2 
jR(3)=rS3 
JR ( 4 ) sr-,24 
JX<l)s.“Xl 
JX( 2 )='-•/? 

JX ( 3 ) s- ; X -3 
J X ( 4 > s r X 4 
O' 15 1 = 1, -.l 
DC 15 J s i , : i X 
l i( I , J) =3. IE- 12 
U X ( I , J ) = ' . i f. - 1 2 
R ( I , J ) = •'. , 1 f. - 1 ? 

T { J , J ) r T I 
T ;: ( I , JJsTI 

n n ( J ) = F 

o ( i ) = i 


, C P L , A :< 5 , A - L , P C S , ? * L 

= • ,F, ' dt-;/(L3*oes f ) * 


r i 


I 


A K 6 = 


I c I 


, / , 1 - X , ' L P L 
t; T’J/( SECHFToQE 


rol 


“TV/ < SEC *FT*Dc5 £)',/» 1 ■ 
' . F , ’ LB/FT»»3 ! ,/> 


A I 


ROS = 


, r l ’ 


co it 



16 


8 5® 
851 

2 v 
21 

121 


A 7 I = 1 1 9, 9 3 
DO 16 ! =1,-1 
T ( I , 1 ) = 7 p 

T I = 0 , 3S?. 

I ' ! P = 1 

AS = 0,893t>>- ! 6 

AL= . a 5 f - P> <; 6 

Ai8'-?:-’i?vt-^/i2 . 

VTsPLQAT ( M 1 ) » A L / a i / C ~ 

WRITE i5.35e») VI 

FORMAT { * f-AX I MUr-i ALLOWABLE VELOCITY \E.' FT/SEC '/) 

GsV'U s; 122 

I TV LEVEL ' , E , / ) 
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22 

33 

32 

C 

42 


dp i re 

' / V 

. * — ■ 

17 1; 

t J y 

i) 

G 


F Q R H .■ 

•T < 

* A 


■ROx 

IMA TE Gq 

T i = r i 

: > ! ■ T 





00 2 1 

. 1 = 

?. : - 

2 



TUA 

■i ) = 

T ( I 

# 

2 5 


oo i; 

L'i J 

= 2, 

1 



T ( 1 , , 

j ) = 7 

(2, 

J ) 



T < ' = 1 , 

C. 

ll 

T ( i" 

2, 

J) 


CALL 

SET 

( R 

* ;, i 

0 | i V ! 

L. • HP ) 

CALL 

VEL 

(E 

0) 



IF ( 1 

IMF. 

/S 

c. • 

2) 

r. 

TO 30 

00 22 JL 

= 2 « 

'•'2 



CALL 

TSP 

< i 

L i 

•10, 

a$» TV. T, 

CALL 

TL.R 

< i 

t 

U 1 


•A L i r i T j 

COOT ! 

; ,-j ij r 






ij 0 T j 4 Y: 

D v 3 2 J L = 2 i ' f ‘‘ a 

CALL TSH (JLi ML i ? •! « t , RH ! ;) 

CALL TL £ I J L . L # T •» » T , R M I . . ) 

CONTINUE 

IL p -l 

PHASE CHANGE CAi.CULA.7 I 0L3 ' . 

C ■ j 7 \j J : J j N t. 

DC- ?:• 1 = ?,! 2 

IE < T N ( I , J) ,LT,T~) CD TO 70 

IF ( ( R ( I , J ) K‘ , ) , GT , HR ( J ) ) GO TO 70 


61 

62 

70 


92 

C 


A H F 

5 

HR 

( 

J 

) 


SUM 

= 

TL 1 

< 

I 

. j) 

_ J { 

SUM 

2 

su 


+ PC I 

i J) 

SA = 

s 

1 i» 


7, 

. 51 

72 

IF 

( 

3 A 

- 

A 

H F ) 

61 

R < I 

, 

J) 

r 

5 

* 5 



194 


T\< I i J)sTF + SA/aHF*3,*' 

GC TO 70 

P( I » J ) = - r -' ( J 5 / 0 , 5 i 7 + i 
T >' I I » J ) = 3 i * T F ♦ ( 3 A ~ A ■- F ) / , 5 1 7 0- 

continue 

n,j ^2 i = 2 1 i i p 
DO 92 J = ? , L 2 
T ( I , J ) = T : , ( I , j > 

PRINT P E 3 L 1 1. T 0 , ( , , , , , . V 

IF <Ti.LT. API) 0 :• TO 2C* 

A T I = A T I + 1. 2 . .0 
R T I = T I / c, 0 . L 

HP I TP (5,194) PT! 

FORMAT ( » T IL-E = ' » F , i Yj.j .,/) 

00 95 1=1,4 
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94 

95 

60S 

96 

97 
132 


8 

12 


13 


13 


00 9 5 J 

j 

= 1*4 




J = J 8 ( I ) 






K = JXUJ 

) 





WRITE ( 

6 

. 94 ) 

J , '•< 

. T < * , 

J-l) , T ( !( , J ) , T ( K , J ♦ 1 ) 

Fj?M..T 

( 

215, 

3 F 1 2 

.3) 


CCK’TIUL 

C- 





vi HITE { 

A 

, 6'"C 

) ( 

1(1). 

1=2,42) 

r o p “i a T 

(1513 

/ 



UiAi I IF (. 

#0 

,96) 




FOF/'A t 

( 

' IU 

P’JT 

FLA G 

»/> 

R Z A 0 ( 6 

I 

97 ) 

FLAG 



r* a r* \i . ^ 

r J **. • 'm I 

(F) 




if (Fla 


Vgf. 

1 oi - 

) Go 

TC 22 


5 TO? 

r ; j o 

SUP 9: OUT ! HE V ZL < ) 

•j I K t S I 0 ) • V(61»3J> , U2 (6ii 33) »U*<Sl#33># M 0 < 8 1 ) 

COMMON /AL^HA/ U, V 
C Q » A •*> H /eETA/ U ? , J f? 

C 0 U •*! 0 '4 / G A P H A / 

00 13 ■' 2 

L s w 0 ( i\ : ) ~ 1 
CO 13 1=2.1 
IF ( L L T , 3 ) GO 
11=1-1 
J J = ■' ,’ 0 ( .' )-l 
IF '(UJ.LT.l) J 
SisFLOAT{ I I 5/FLG' T< JJ) 

IF (Si,GE,v.,999> 50 TO 12 
CO 3 J = 2..'-2 


Pi, N 2,'s3, Ui ,:-'2,M3 
TO 12 


J '< - J - 1 

S 2 = F L 0 A T ( J K ) / F L 0 . T < v ? ) ' 

IF (S2.LT.0i) GO TG .. 

uz< m, n =u(p, j) 

uau, n = v( - , j) 

GO TO 53 
CO>:T[NL£ 

GO TO 13 

JZ< -!. I )=, V .1L-12 

uk( D=n.iF-i2 

COi-'TIHlit 

retofl 

E H 

SUE- ROUTINE VIVV( '-MAX) 

CIMT- S!0.\ 0(61.3/) >”(51, 30) 
c 0 M M 0 U / A l F -i A / u . V 
C 0 * '* 0 N / C A M H A / •'•! T , N 2.^3,: 

2 1 = F L 0 A T ( ~ 5 2 ) 

VVs-1.0 
8 = 1,? 

A*-»i. 2/'SGRT(3.? ) 

DC 10 1=2. h 2 

Y V = V •/ + P , G i P X 
C = Y -/ / 9 

9A= ( C * » 3 . TfC)/ ( a»*3. >A) 

VV( I ) = >• : ' A V A 
? I = 3 , 1 a 1 5 9 6 2 

Y = •/ . 

00 15 I = 2, HP 

Y = Y + J. .'C-./Pi 


1 


1 . 


, V V ( 3 ) , SF(3'/i) ,03(30) 
2 » f -* 3 


V8(30) . VO (33) 
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15 

18 


32 

31 


36 

40 


53 
6 3 


X c \ , Z 

sr ( i isct-wt? 

UR ( I. ) s^Sl -'-’H ( P I*X > C S < P ! *Y 5 «P I 
VS ( : )sC0Sri< P ! »X)*3I.V<PI«Y) *p I 
X = 2.0 

DO 1 8 ~ J = 2 » 2 
X = X ♦ 1 , f/* ./ Q i 

VO ( «l + l - J ) = " I *CC3H <PI*k)«S!-‘<i(Pl»Y) 

YsC.5y.-l.'C/*l 

Jl=('-4 + U/H 

00 36 M = -J 1 » > ". 2 * 1 

Y s Y + 1 , r / R i . 


y y 0 0 1 = C # 1 ■ £ 

X = X + 1 V 0 / Pi 

SFT = SUH(Pt#X)«StN<?!«Y) 

Uls-P I »s I ( ? I #X ) «C05 ( P I *Y ) 

Ll»M?/2 
00 3;. !=1,L1 
K X s • < 5 1 - 1 

IF (SFT,GE,SF<KK)> GO TO 33 

GO TO 31 

CONTINUE 

K K = J l 

U H A X = u 6 t K K 1 
U BO- V V ( X ) 

U < !' , M) =Ul«JL ,n /i.'iX \x 

V;.sP I acO^'i < ~ I*X) I r l** ) 

V ( iV v ) s ( 1 , » - .A S 5 ( V / V "i ( KK ) 5 ) #U 8 0 

IF ( A Q c ( v / • / R ( vt * ) ) , OF. ,1.02) V ( M , M ) :2. 3£1*U&0 
continue 

0 0 4 3 .'••• = 3 I w 1 
00 4 /■ ‘2 

u ( ••'i # ‘ ) S-U < Yl-'M 

V ( IV , " ) = '••• < <i , 1- P ) 

COOT I UUP 

f-Ml s. 2 
MR2 S ,h1-1 
MP3='P?-1 
00 5-2 i v = 2 1 2 

DO 5 ■ Msi-t j, » MRi 

V (0, M )■•-}, lE-12 

U { K j . •! ) s V V ( 0 ) 

IF (K.e).Ul) U(X,M) =0,1 £-12 

cortimue 

00 60 f- = 2 , M 2 
DC 8- i'-=-Rl,03 
MMsM-MP 2 

U ( P 1 — : v » ; l ) ~ u ( ■>< - N i '•< .> i ” ) 

V ( ‘oi-\r. , H ) =-v ( M — P.3 , • . } 

C 0 ' : T i N I.! £ 

RETURN. 

C ‘-r; 


S U 5 P 0 U T 1 : 1 E SET ( o , 0 . H L , - 1 •< J 


0 l r 

••e v s 

i or- max,- 

5 1 1 » '•) ( a i ) . ■ 

JL<3GI 

COr 

•-•ION 

/GAiU-Pt/ ;■ 


M2. M3 

00 

1 ? 

I=2.P2 



0 0 

6 J 

*2.*: 2 



IF 

(P ( 

! . J)*0,5l7'.GE.-iR< j) ) 

GO TO 



8 

12 ! 


18 

20 


10 


12 


15 

20 


N0< I )sj 
GO TO 10 

continue 

MO ( I ) s f ' 2 
CONTI >\VZ 
DO 22 J«2#n2 
«L( J)s2 
DO 18 I - 2 » : ! 2 

IF ( R C I » J>*0.5i7.'LT,»:R( J) ) GO TO 10 
HL ( J > = I 

GO TO 20 

c o m r i N i* e 

MLf 

c 0 f- ’ T I r.jU£ 

RETURN] 

Z>*C 
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S U n 

R 

P 

u 

TU 

r* 

T 

SR { 

? , 

NO 

l 

A 

A 

s, 

T 

0 ! M 

£ 


Slow 

T 

s : 

< 81, 

3 2 

) . 


T 

( fi 


OIK 

r- 

c. 

jv 

r* 

O 

!OM 

u 

( 

3 ?) , 


3 

J 

» 

Y ( 

•J 

COM 

H 

0 

* j 

i j 

/si/ 


a , C , 

» 

i.i 

# 

V 




COM 


0 

* i 
! x 

/■> 

2/ 


D ’3 i D 

« 

Or 





COM 

^■j 

»./ 

V 1 

/G 

A • 1 

M 

A / N 

1 # 

'■■>2 

1 

V 

3 , 

;-t 

IF 

( 

V 

0 

( I ) 

.0 


.M 2 ) 

r 

0 

7 * 

0 

1 

-I 

IF 

( 

*. 

n 

( I ) 



* / 2 5 

G 

0 

T 

r*1 

2 


! p 

< 

* , 

n 

( .1 ) 

. £ 


, M 3 ) 

5 .J 

0 

r 

i 

n 

2 

y 

O' 

(1 

y 

..i 

{ 

I 

) 










R A o 

= 

r 

L 

M A T ( 

- 

1 ) » i." 

+ 

R M 

I 

is’ 

/l 

2 

A ( K 

\ 

/ 

= 

m 


:*j T 

/ 

o a / n 

4 

'* f 



A 

is 

0*1 

« 

j 

+ 

c % 

« A 

- J 

* -j t / 

r ? *y 

/.j 

s 




C(K 

) 

= 

■w 

A T* * 

QT/ 

0 R / D 

r • 

•• • 

'/*, 

3 


A 3 

n 

0 < K 

) 


T < I , 

/ ) 

♦ 

A 3 * r- 

T/ 

r\ 7 

v it; 

/ 

Q 

H « 


0 ( K 

) 

- 

~ I 

( ) 

•• r 

7 

K ) o T 

( T 

• K 

- 

-L 

) 



K K S ^ + 1 

DO Vi J = K K ( ,'i -3 

R A 0 ~ F L C 4 T ( J - 1 ) * 0 ~ ♦ R " : I i / 1 2 • ■' 

M J ) s - L > * 0 T /0 s !/ ,0 ♦ 0 i ' * R S # 0 T / :T A 0 / 0 9 
C ( 0 ! -»• j C'^v 1 / •■t / : J ' fl *7 O A # j T / ■-*. L f; / ;*) h‘ 

0 < J ) = T ( I , J ) + A S * 0 T / D S / 0 2 » ( T ( I + 1 * J ) ♦ T ( I - 1 » J ) - 2 , 0 * T < I > J ) ) 

JsM2 

R A D s " L r ' i T ( J- 1 . ) -OR + R:- T h/'V'. , 0 

a ( j ) s » a j # j t / u r / o i ♦ 0 , o # a s»riT /R ad / o R 

J T / “ U iJ / i R 

0 ( J ) s T ( I , J ) + A 3 » 0 r / fj 2 / rj -2 # ( t ( j * 1 , j ) + t ( I 1 * J ) - 2 , * T < I , J ) ) 

0 ( J ) =D v J) -C ( J.) *U l , j*D 
W ( s< ) = 9 

Y ( K 5 s f: ( K ) / » { K ) 
no 12 J s K K , M 2 
5.5 s C ( J-1 ) /•■• < J-1 ) 
w(J)s3*A(J)»98 

Y< J) = (P! j)-a( J > * r < J-1))/-. ( J) 

T\ < I , y.2 ) = Y( ; 2 ) 

DO 15 JSK3.K.-1 
L = J 

80 = 0 ( L ) /’■•• ( L ) 

T \ U , L ) s Y ( L ) -QO*T.j( I , L*1 ) 

GO Tj i j'? 

R A psFLC a T ( M-i ) s n R -r R H I •' ■ / 1 2 , 7. 

A a = - A S ^ 0 T / D ■' / J >’• ♦ •: , 3»GT°AS/:R a P/PR 
8 B 3 i • S * 2 . /,■ / S » 0 t / n R / ; « 
c C = A A ^ A s * u T /RAM/ f j R 
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25 


13 ! 


13 


OD - T ( 1,m?)+AS»0T/D2/02*<T ( 1 *1 » "'2 ) ♦? { ! -i , !V i? ) -2 , *T ( I i M2 ) ) . 
T ■■< [ , :\2 ) s(00-AA*t< ! » i ; 3 ) - 0 C * T ( I , M i ) ) ✓ -3 S 

GO Tv 100 

R A 0 s T L 0 a T ( -3 > • GR + *** ! ‘Vi 2. 3 
R a = - AS * 0 T / 0 :? /OR* • , 5*aS* 0T / ?i A 0/ D R 
C C = R A - A S & 0 T / Q R / R 0 
R A [) ~ P a r. - 0 r? 

RC = -AS«OT/UR/OR- ,5*/<3»0T/RaD/0R 
A A = R C ♦ A S » 0 T / R A 0 / : • R 
6 8 = 1 , 3 * 2 - ■- * A S * f < T ./ 0 R / C R 

01 = T £ 1 , ••3>*A$»f’T/oZ/e2«(T< 1 * 1 , H 3 ) * T ( I - 1 . 03 ) -2 , «T ( 

Olrul-A a*T f J , ' 3 - ) 

0? = T( !,•:>:) ♦aS*OT/dH j /•”.**( T( I*1 ,m2)+T i I - I,:- -2 ) -2 . »T ( I , M 2 ) ) 

02sD2«CC*T < ! » 01) 

rsni/BP-Q2/ r -A 

G = f!c/Qc-U"/RA 
T.’*{ I f n2)sr/G 

T'v ( ! ,MS ) =Qi/83-R0*T'V < J ,m-m/33 
R £ T ! j R M 
E 0 

SUBROUTINE TIR < I , NO , At , IN , T , kk ! M > 

0 I t NS I 0M T a ( a i , , T ( 6 j. , 3 3 ) # N C ( 3 1 ) . A < 3 2 ) # C ( 3 0 ) » 0 2 < fi 1 » 3 0 ) 

01 MENS I ON UR (3 1,30) 


0 I M 

£!• SI 

*. A 
V *t 

0(32), 


(30) , Y(32) 

COM 

M pi f» i 

/s 

1 / A , C , D 

• 

v, Y 

C'j;' - . 

MOO 

/s 

2/ DZ,D 


I*- 1 ! • 

CO- 

: i o n 

/p 

ETA/ 02 

# 

: • P 

M 

COM 

••nn; 

/d 

A MM A/ 


, i' £1 , :'!C , J. , 

L>J = 

.'0(1 

) - 

J. 



! r 

( L M . 

LE 

| 3. « * * 4 ) 

r 

>1 i a 

IF 

(Lu'i 

EO 

• ?. ) C> 0 


0 2 0. 

IF 

(UJ'i 

r r> 

,o) n- 

T 

r: p o 

n A 0 

2 fA y\ J 

;:/i2,2*pr 






A< 2 )si*:R( I » 2 ) »:-T/ 2 t 3 /O^*'AL* r )T/r.R/OH* f l« 5 *AU* 0 T/RAO/ 0 R 


1 < < 


1 * a i « r, t / 3 v n r 


C ( 2) S3, S*u* (1,2) oijT/OR- AL^OT/OR/OR-.?, . 5* a L *{.'T /RAO/OR 
0(2 ) = .*, 3» 3 7*112 ( ! . 2 ' • ( T< t- 1 , 2 ) -T < !* 1 , 2 > )/ -Jr + T ( ! , 2 ) 

0 ( 2 ) s.0 ( 2 ) ♦ Ai,*DT* ( T U +1 , 2 ) ♦ T ( I -1 , 2 J -2 ,'^*T < 1 , 2 ) ) / OZ/OZ-A ( 2 ) *T U 1 1 ) 
L C = L ;v - 1 
0 j l J - v'. • L- C 

RaOsFLCAT ( J-l ) «0:S*R- v ! N / 1 2 ./> 

A ( J )=-•;• , 5«UR( I, j ) « 0 T / D « - A L » n T / 0 R / 0 R ♦ 0 , 5»AL*0T /RaO/OR 
C ( J) =0 . ?»!;- ( I , J) » 0 T / •“■ r < » a u » ; i T / D R / 0 R - ^ , 5 <> -A L * 0 T / R A 0 / 0 R 
0 i j) = 2. S*OT«UZ ( t , j) » ( T< !-i. J)-T( 1*1, J) )/r.-H*T ( I , J) 
r, ( J)sD< J ) + AL* r )T*<T( 1*1, J) + T« 1-1, J)-2.0*r< I , J) )/DH/OH 

RAOsTlOaKU !)*^R + RMr'/12 t •' 

A ( L 'i ) s • 7 , 3 ^ ‘ ^ ( I , L ' i ) a 0 I / v) i - A L 0 0 T / 0 R / fj R + •? , *■ *• A L *0 T / r! A 0 / OR 
C ( L '•!) = :•, 5 » *-■ ( lit, - ) °0'/ y R - A U * 0 T / 0>- / 0 R - a , ~ * a l a 0 T / R A 0 / D R 
0 ( U •'* ) * 0 • 5 *0 T * ‘j 3 ( • , L • » ) » ( T < t - 3 f l, f.: ) - T < I * 1 , L ' • ) ) / 0 H + T ( I , L M ) 


0 ( L ■' -i 

) S n 

(Li) 

* A L 

C(L*i 

) = 0 

(L.M) 

~C ( 

( 2 ) 

•• :~\ 
« r X 



¥(2) 

= 0 ( 

2 ) / v 

(2 5 

no ,i 

2 J 

s3, L 


50 sC 

( J- 

t ) / 

( j- 

!'• ( J ) 

= B- 

A (J) 

*35 

Y ( J ) 

= <!? 

( J>- 

A ( j 

r. ( ! 

• L ; - 

)=Y( 

L ) 

DO 1 

5 J 

=10, 

r>. 

C l - 


L = 


12 



BB = C«U/v:(L) 

15 T\< !,L) = V!L!-iie«T-;( I , L + l 3 

GO TO 1 ?.?. 

20 RACsSHI VlS.S+CR 

» A 5 - A L * 0 T / c ? / 0 F- ♦ ,5*..U*DT/>"A0/DR 
P8*l « 2*2, 0»Al>'~T/DR/ r -» 
f< C = h A « A L * 0 T / > A A / ■ ■ 


25 


10- 


5 


*0 

= T ( I , 2 ) * A 

i « 0 T / 

•:•=/. o z 

o ( 

T< 1 + 

1,2) + T( I 

-1.2) 

-2 . *T ( 

I 

.2)3 

TN 

( 1,2): 

5 ( * f ‘ 

-PA«T 

Mil) 

— > K 

C ft T ( 

! .3) )/ 

Re 





GO 

TO lv 












n u 

QskM ! 

via 

.0*2. 

7 » 0 A 









R A 

- _ . ' ‘ * r. , 

*• • - • 

HjRC 

I » >-*' •) 

»DT/C- 

’ A *• 

A L » 0 

T/QR/J 

h + 

3 . 5 » A L * 0 T / R 

AD/DR 

RC 

c >i a + i; : 

<( I , 

!j -: ) *n 

r / 0 0 *. 

A L 

» r • T / 

RaO/QB 






R3 

= 1 . 0 * 2 V 0 * 

A L ft 0 T 

/DR/O 

R 








2 A 

r» - ^ a [* * 

. . r ^ 
■ i : r( 











Cr- 

s3 , 5«>. 

jq ( i 

i L *“ 1, 

) oOt/ 

; 1 CO 

- A L * 

ot/pr' 

Or< 

-a. 5* 

AUftOT/ 

R 

A D/OR 

ab 

= - ’ 7; , r - # U n ( 

I . L-- 

) * 0 T 

/‘D 

F-AL 

ft 0 T / Q y 

/O 

k + 

ft 0 T « A L / 

RaO/OR 

01 

= 0 , ? « 

I 

. LO-i 

j «• 0 T / 

02 

»<T< 

! - 1 , L ‘ 

- 1 

5 - T ( I 

+ 1 .LO- 

1 ) ) + T < I ,1-4-1) 

01 

= 01+ t>l »DT 

*(T(1 

+ 1 . L : ■ 

4 

"* Jl 

> + T< 

t-l.L:. 

4 

) -2 , ft 

TH, L 

- 
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